% 使用swt计算小波边缘提取
% P 要提取的图像
% scale 最大尺度
% threshold 阈值，常数，需要修改
% wavelets_index 使用哪种小波
% E 返回的边缘提取后的边缘图像
% M 返回的边缘提取后的模值图像
function [E M] = compute_wavelets( P, scale, threshold, wavelets_index )

% 为了与conv版本的scale保持同步
scale = scale+1;

%% 将图像的长宽置为 2^scale 的整数倍
[C,L]=size(P);

if mod( C, 2^scale ) ~= 0
    C = C - mod( C, 2^scale );
end
if mod( L, 2^scale ) ~= 0
    L = L - mod( L, 2^scale );
end

B = P(1:C,1:L);
P=double(B);


%% 数据定义
% 幅角图像
Arg = zeros( C,L,scale );
% 模值图像
Mod = zeros( C,L,scale );
% 求出局部极大值后的图像
Mod_after_max = zeros( C,L,scale );
% 将Mod_after_max中的非零值置为255
Mod_max = zeros( C,L,scale );

if wavelets_index == 1
	% 正交二次样条
	h = [0.125,0.375,0.375,0.125]; 
	g = [-0.9324,0.6736,0.2253,0.0335];
	[A,H,V,D] = swt2(P,scale, h, g);
elseif wavelets_index == 2
	% 正交三次样条
	h = [0, 0, 0.0625, 0.2500, 0.3750, 0.2500, 0.0625, 0, 0, 0];
	g = [-0.00008, -0.01643, -0.10872, -0.59261, 0, 0.59261, 0.10872, 0.01643, 0.00008, 0];
	[A,H,V,D] = swt2(P,scale, h, g);
elseif wavelets_index == 3
	% 非正交二次样条
	h = [0.125,0.375,0.375,0.125]; 
	g = [-1.0,1.0]; 
	[A,H,V,D] = swt2(P,scale, h, g);
end


%% 计算幅角图像和模值图像
for i=1:scale
	Arg(:,:,i) = atan( V(:,:,i)./H(:,:,i) );
	Mod(:,:,i) = sqrt( V(:,:,i).^2+H(:,:,i).^2 );
end

%    -------------
%    | 1 | 0 | 3 |
%    -------------
%    | 2 | X | 2 |
%    -------------
%    | 3 | 0 | 1 |
%    -------------

%% 计算局部最大值，并进行非极大值抑制
for i=1:scale
	Mod_after_max(:,:,i) = do_modmax( Mod(:,:,i), Arg(:,:,i) );
	for x=2:L-1
		for y=2:C-1
			if Mod_after_max(y,x,i) > threshold
				Mod_max(y,x,i) = 255;
			end 
		end
	end
end

%% 连接边缘
for i=1:scale
	Mod_max(:,:,i) = do_edgeconnect( Mod_max(:,:,i), Arg(:,:,i) );
end



maxvalue = max(max(Mod(:,:,scale)));
Mod(:,:,scale) = Mod(:,:,scale) * 255 / maxvalue;

M = uint8(Mod(:,:,scale) );
E = uint8(Mod_max(:,:,scale) );

end

